Transport Phenomena in Fluids: Finite-size scaling for critical behavior 
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Results for transport properties, in conjunction with phase behavior and thermodynamics, are 
presented at the criticality of a binary Lennard- Jones fluid from Monte Carlo and molecular dy- 
namics simulations. Evidence for much stronger finite-size effects in dynamics compared to statics 
has been demonstrated. Results for bulk viscosity are the first in the literature that quantifies crit- 
ical divergence via appropriate finite-size scaling analysis. Our results are in accordance with the 
predictions of mode-coupling and dynamic renormalization group theoretical calculations. 
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Understanding the properties of fluids is important 
from both basic research as well as technological point 
of view. Particularly, fluid behavior in the vicinity of 
critical point poses many interesting questions of fun- 
damental importance [l|-l8| . Fluids with short-range in- 
teractions, exhibiting gas-liquid and liquid-liquid tran- 
sitions, are expected to have the static critical expo- 
nents p = 0.325, 7 = 1.239, v = 0.63, a = 0.11, re- 
spectively, for order-parameter, susceptibility (x), cor- 
relation length (£) and specific heat, thus belonging to 
the three-dimensional Ising universality class [7[. On 
the other hand, it is expected that model H 0] should 
define the dynamic universality class for both the tran- 
sitions. In dynamics the quantities of interest are shear 
(rj) and bulk (() viscosities, thermal diffusivity and its 
analog, the mutual diffusivity (Dab) in a binary fluid, 
the critical singularities for which are given by [2|-|4( 
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e(= \T — T c \/T c ) being a measure of the temperature 
(T) deviation from the critical value (T c ). In Eq. JT]), 
exponents obey the scaling relations [2| 



Xd = 1 + Xr), 



x n = 0.068, z = 3.068, (2) 



where z characterizes the divergence of relaxation time 
t, at T c , with system size L as 



(3) 



In addition to other static and dynamic properties, 
particular focus of this work is to understand the critical 
behavior of bulk viscosity that describes the response 
of a fluid to a compression or expansion. Even though 
the study of bulk viscosity is thought to be important 
for compressible fluids, the theoretical prediction, as 
our simulation results will also reveal, for similar criti- 
cal enhancement for both gas-liquid (compressible) and 
liquid-liquid (incompressible) transitions is certainly in- 
teresting. While there are experiments 0] probing the 



critical behavior of C, simulations are rare de- 
spite this being very important in the description of 
the damping of longitudinal sound waves. Dyer et al. 
13| in fact pointed out the difficulty of studying bulk 
viscosity and suggested the need of significant effort to 
understand dynamics of continuous model fluids. The 
only noteworthy study of bulk viscosity in the context 
of criticality, so far, is due to Meier et al. [ll| for gas- 
liquid transition of a single component Lennard- Jones 
(LJ) fluid who, however, did not quantify the critical 
divergence. While their data suffered from large error 
close to T c , their observation of strong enhancement of ( 
far above T c (~ 4.5T C ), as is also observed by us, due to 
extremely slowly decaying pressure fluctuations, is very 
interesting. In fact, to the best of our knowledge, ours 
is the first simulation study of bulk viscosity that quan- 
tifies its critical divergence. In addition to confirming 
the theoretical predictions for critical exponents 0-0] > 
this work also provides direct evidence for stronger size 
effect in dynamics compared to statics. 

Apart from understanding the universality, computer 
simulations have been instrumental in providing many 
other details as far as static properties are concerned. In 
contrast, simulations of critical dynamics are very rare 
\.5\ the primary reason for which being the criti- 
cal slowing down, as embodied in Eq. ((3]), that brings in 
additional complexity for the computation of dynamics 
over statics where finite-size effects are the only diffi- 
culty. Also, for molecular dynamics (MD) simulations 
in micro canonical ensemble, which is needed for perfect 
preservation of hydrodynamics, it is extremely difficult 
to control the temperature to the desired value for a 
prolonged period of time due to truncation error. This 
may be a necessity at temperatures close to T c because 
of the presence of long-time tails 0, HH, particularly 
for quantities showing strong enhancement. All these 
problems combined together is suggestive of avoiding 
brute force method of simulating larger systems close 
to the critical point. The finite-size scaling method em- 
ployed here to understand the results will demonstrate 
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that if appropriate strategy is devised, all these hurdles 
could easily be overcome. Apart from critical phenom- 
ena, such methods could be useful in the study of other 
slow dynamics, e.g., glassy dynamics, where tools from 
critical phenomena are being recently adopted to un- 
derstand growing dynamic length in supercooled liquids 



We use a binary fluid (A+B) model [1J,[18| where, for 



l\) < r Cl particles at fi and fj interact via [1 
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U{r c ) - (rtj - r c ) (dU/drij) ri _ r 
0, with U{r l0 ) = 4e Q(3 [(cr/r l /) 12 



(&/ r ij) ] i a - ft G A,B) being the standard Lennard- 
Jones (LJ) potential. For the choice Eaa = Ebb = 
%£ab = £, we have a fully symmetric model that 
gives a demixing transition at 14] T* — ksT c /e = 
1.4230 ± 0.0005, for r c — 2.5a. Phase diagram of such 
a system can be obtained from a semi-grandcanonical 
Monte Carlo (SGMC) simulation [13], where in addi- 
tion to standard displacement trials, one introduces 
identity switch (A — > B — > A) moves, thus allowing 
for fluctuations in concentration x a (— N a /Y2p Np) of 
species a. The distribution P{x a ) of concentration 
fluctuation has double and single peak structures re- 
spectively at temperatures below and above T c . While 
from the location of the peaks below T c one can obtain 
the phase diagram in xa — T plane, static concentra- 
tion susceptibility (\) above T c can be calculated as 
k B T X = X*T* = N{(x a 2 ) - 1/4), where the term 1/4 
corresponds to a critical concentration x c A = 1/2 dic- 
tated by the symmetry of the model. 

Transport properties were studied via MD simulations 
in microcanonical ensemble with £ and the Onsager 
coefficient Jz?(= X^ab) being calculated from Green- 
Kubo (GK) relations jl& [U 
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dt(J™(0)J™{t)}, 
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where a xx = a xx - P, a xx (= J2iLi[ m i v ix v ix + 

\ Y^j (%i — x j)Fxj\) being the diagonal elements of the 
stress tensor with P = (cr xx ) (GK formula for 7/ in Eq. 
(|4| contains the off-diagonal elements of stress tensor) 



and J, AB (t)(- XB^iAviAt) ~ ^aES^bW) is thc 
concentration current with v^ a (t) being the velocity of 
particle i of species a at time t. In Eqs. (HJ) and (0, 
V is the volume of simulation box, m is the mass of a 
particle and to[= (ma 2 /e) 1 ^ 2 ] is the LJ time unit, which 
we set to unity. All results for dynamics are obtained 
from MD runs, with integration time step At = 0.005, 
in a periodic cubic box of length L, in units of <r, after 
averaging over 160 independent initial configurations. 
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In Fig. QJa) we show the results for x as a function of e 
for L* — L/a = 10 and 18.6. The continuous line there 
has a power-law form with exponent 7 being fixed to its 
Ising value 1.239. While this confirms the Ising-like be- 
havior, the consistency of the data for L* = 18.6 with 
the solid line over the whole region is suggestive that 
finite-size effects did not appear yet. In view of the fact 
that finite-size effects were pointed out to be stronger 
in dynamics and a very strong background contribu- 
tion was found in the study of mutual diffusion [14| , we 
revisit it in the following. 
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FIG. 1: (a) Log-log plot of x vs. e, for L* = 10, and 18.6. 
The solid line represents the critical divergence of x with 
exponent 7 = 1.239. (b) Plot of Onsager coefficient AJt?/T* 
vs. e for L* — 10, and 18.6, on a log scale. Ths solid line here 
has critical exponent v\ = 0.567 and amplitude Q = 0.0028. 

In Fig. Wijo) we study the critical enhancement 
AJS?(T)[= % - JS? 6 ] of Onsager coefficient, with Jgf 6 be- 
ing the contribution coming from short-range fluctua- 
tions and needs to be taken care of far above T c , where 
critical enhancement is small. In the following we treat 
Jz?b as a constant, albeit weak temperature dependence 
that it might have. Here we plot A££ (T) which has the 
expected critical divergence 



AJf = QTe- 



v x = 0.567, 



(0) 



as a function of e, by adopting the constant value of 
Jz?b = 0.0033 as obtained in an earlier study [L4|. Upon 
imposing [l4| Q = 0.0028 and v\ = 0.567, good consis- 
tency of the solid line is obtained with the simulation 
data, for large e. This, in addition to directly confirm- 
ing the theoretical predictions as well as the conclusion 
drawn from the previous finite-size scaling study [LJ], 
with very limited data, regarding the exponent and am- 
plitude, is also indicative of a rather wide critical range. 
On the other hand, it is interesting to note from the 
comparison between Fig. (Ha) and Fig. Gib) that s i ze 
effects are appearing much earlier in Jz? than in x, which 
requires appropriate attention to understand. With the 
knowledge about the spread of the critical region, we 
move forward to devise a strategy for a finite-size scal- 
ing analysis [a, LL9j, |2l| , that will require only small sys- 
tems and large temperatures so that difficulty due to 
long-time tails could be avoided. This will be tested 
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FIG. 2: (a) Phase diagram of the model in xa — T plane for 
two system sizes. The filled symbol corresponds to Tj? for 
L* = 8, while the cross (x) locates T c °° = T c . The dashed 
line there is a fit to the form m = \xa — 1/2| ~ e' 3 , taking 
data close to the critical point and unaffected by finite size, 
(b) Plot of vs. 1/L* where the continuous line is a fit to 
the form (0 with v = 0.63, including only the four largest 
system sizes. 



with the better understood quantities, \ an d -Sf first, 
before applying it to f . 

As a first step, in Fig. [2fa) we show the phase behav- 
ior of the present model for different values of L that 
exhibit strong size effect close to the critical point. We 
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as the tem- 



define a finite-size critical point 
perature where P(xa) in the SGMC simulation gets a 
single-peak structure from a double-peak one with the 
increase of temperature. This is represented by filled 
symbol for L* — 8. Note that true meaning of a crit- 
ical temperature can be assigned only when L — > oo, 
which for the present case is marked by a cross and was 
obtained [l4| in an unbiased manner from the method 
of intersection of Binder parameter 23] . In Fig. &h) 
we demonstrate the variation of T C L with L. The con- 
tinuous line there is a fit (including data only for four 
largest values of L*) to the expected scaling form in the 
large L limit, 



(T - Tj 
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(7) 



The deviation of data for L* = 10 and 8 are due to 
corrections to scaling for small values of L, which can 
be numerically accounted for [!, Q by replacing L* by 
L* — I* in the abscissa. Indeed for I* — 2 we get a 
perfect fit passing through all data points, which, in 
fact, has been used to obtain for intermediate L 
values. Nevertheless, even for L* — 8, correction is 
very small and the simulation data deviates from the 
solid line by less than 1% which is negligible compared 
to the thermal fluctuations during the MD runs. 

At this stage, we define an effective finite-size critical 
point [H from T C L [= T C L (1)} as 



T^(f)=T c + f(T c L -T c ), 



(8) 



which has the same power-law convergence to T c as ([7| • 
One can study critical behavior along different /-loci as 



a function of L, when, for an observable O (~ e x °), one 
obtains the scaling law 
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where the amplitude will depend upon the value of /. 
Fig. [3] demonstrates this for \ an d Aj£? for two values 
of /, where we have plotted \ vh an d {^Y /vx vs L. A 
linear behavior upon imposing v — 0.63, 7 = 1.239 and 
v\ = 0.567 validates this strategy. Note that largest 
value of / considered here is 65 which gives T C L (65) = 
3.875 for L* = 8. Nice consistency of the data for whole 
range of L is suggestive of only weak corrections even 
for the smallest L considered. 
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FIG. 3: Plots of x vh and {^Y /vx vs. L* along different 
/-loci. Data for different quantities and / values have been 
appropriately scaled to collapse. The continuous straight 
line is a guide to the eyes. 



Having demonstrated the usefulness of such method, 
encapsulated in Eqs. ([8j and (j9)), we adopt it to quan- 
tify the critical divergence of £. Due to the technical dif- 
ficulties to calculate it at lower temperatures for larger 
systems, in Fig. Ufa) we present it only for / = 65. Very 
linear look of the whole data set on a log-log plot is sug- 
gestive of only small background contribution. Signifi- 
cant increase of £ over only small range of L* £ [8, 12], 
signals a strong divergence. The continuous line there 
is a fit to a power-law form ~ L x< giving ~ 2.96 
which differs only by 2% from the theoretical prediction 
(almost indistinguishable dashed line) z — ^ ~ 2.89. 
Note that more recently Q it has been pointed out 
that the exponent is closer to z. While this confirms 
the expected theoretical behavior, we estimate the non- 
universal critical amplitude from the following exercise 
which will also provide a more direct confirmation of 
the exponent. In Fig. QJb), we plot £ as a function 
of e. Here, from our experience with Jzf, we choose a 
range with e > l(that includes the last four points) for 



a fitting to the form £ 



by fixing xq to 2.89, 



which gives a critical amplitude Aq = 6.6 ± 1.0. Here 
the point of deviation of the simulation data from the 
solid line is consistent with the appearance of finite-size 
effect in Jzf. While the solid line provides an excel- 
lent fit to the selected region, an effective, though much 
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smaller and misleading, exponent could also have been 
obtained from a fitting to the whole data set which has 
an average linear look on log-scale. 




FIG. 4: (a) Log- log plot of ( as a function of L at T C L (/) 
with / = 65. The continuous line is a fit to ~ L x z giving 
X( = 2.96 while the dashed line corresponds to an exponent 
2.89. (b) Plot of C vs e for L* = 10. The solid line is a fit 
to the form A ( e~ 182 giving A c = 6.6 ± 1.0. 

In summary, dynamic critical phenomena is studied 
in a symmetric binary fluid. Consistency with pre- 
dictions of dynamic renormalization group and mode- 
coupling theories has been established. Quantitative 
understanding of the bulk viscosity via computer sim- 
ulation is the first in the literature. Critical region ap- 
pears to be rather wide so that with appropriate appli- 
cation of finite-size scaling method it has been possible 
to stick to only small systems at large temperatures. 
A possible reason for stronger finite-size effects in dy- 
namics compared to statics could be back-flow due to 
periodic boundary conditions — however, significant at- 
tention is required to settle this important issue. 
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